Gutzwiller approach to the Bose-Hubbard model with random local impurities 
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Recently it has been suggested that fermions whose hopping amplitude is quenched to extremely 
low values provide a convenient source of local disorder for lattice bosonic systems realized in 
current experiment on ultracold atoms. Here we investigate the phase diagram of such systems, 
which provide the experimental realization of a Bose-Hubbard model whose local potentials are 
randomly extracted from a binary distribution. Adopting a site-dependent Gutzwiller description 
of the state of the system, we address one- and two-dimensional lattices and obtain results agreeing 
with previous findings, as far as the compressibility of the system is concerned. We discuss the 
expected peaks in the experimental excitation spectrum of the system, related to the incompressible 
phases, and the superfluid character of the partially compressible phases characterizing the phase 
diagram of systems with binary disorder. In our investigation we make use of several analytical 
results whose derivation is described in the appendices, and whose validity is not limited to the 
system under concern. 



I. INTRODUCTION 

Since the seminal paper by Fisher et al. [1], disordered 
bosonic lattice systems have been the subject of active 
investigation. The recent impressive advances in cold 
atom trapping allowed the experimental realization of the 
prototypal bosonic lattice model, i.e. the Bose-Hubbard 
model [3|- Different techniques have been devised for 
the introduction of disorder in the system |3|, such as 
speckle field patterns 4], incommensurate bichromatic 
optical lattices l5] , localized fermionic impurities ijj] . In 
particular Ref. Q provides experimental evidences of the 
hallmark phase of the disordered Bose-Hubbard model, 
i.e. the compressible and non superfluid Bose-glass ili]. 

At the theoretical level very diverse techniques have 
been employed in the study of the disordered Bose- 
Hubbard model. A non exhaustive list includes field- 
theoretical techniques [l], 01 quantum Monte Carlo simu- 
lations 3, mean-field schemes 0, [13, [H E [IS H IH 
and others [ll,[i3- 

Here we are interested in the case of fermionic impu- 
rities. Bose-Fermi systems have been studied by sev- 
eral Authors, ^M,- If the kinetic energy of the fermionic 
atoms is negligible, e.g. due to a strong suppression of 
the relevant hopping amplitude, the impurities localize 
at random sites of the optical lattice [a, [l^. The sys- 
tem can be hence described by a Bose-Hubbard model 
with random local potential characterized by a binary 
distribution. Several years after an early discussion of 
this model [20] the features, new features induced in the 
phase diagram of the system by binary disorder has been 
discussed in recent Refs. [2l|, [23, l2j, [2J| . A characteris- 
tic feature of such a phase diagram consists in the pres- 
ence of noninteger-fiUing incompressible lobes. Mering 
and Fleischhauer [22] provide simple arguments showing 
that the phase diagram of the disordered model does not 
depend on the impurity density and can be straightfor- 
wardly derived by that of the homogeneous case, at least 



as far as compressibility is concerned. In particular one 
can recognize fully compressible, fully incompressible and 
partially com,pressihle regions. While the first and the 
second are clearly superfiuid and insulating, respectively, 
the question arises about the superfluidity of the par- 
tially compressible regions, at least on high dimensional 
lattices. Indeed, as discussed in Refs. [22, [2^, the par- 
tially compressible phase is bound to be insulating, and 
hence Bose-glass, on ID lattices. 

In this paper we describe the zero-temperature mean- 
field phase diagram of the Bose-Hubbard model with 
binary-distributed disorder. First of all we show that the 
above compressibility scenario independent of the impu- 
rity density [24] is confirmed also by our site-dependent 
Gutzwiller approach. Moreover the analytical tractabil- 
ity and the computational affordability of this technique 
allows us to investigate the superfluidity of the partially 
compressible phase both in one and two dimensional sys- 
tems. In particular, we address the issue of quantum 
percolation which, as already pointed out in Ref. [6|, is 
expected to play a crucial role in this problem. While we 
confirm that on one-dimensional systems the partially 
compressible phase is substantially insulating, we find 
that in higher dimensions the system always exhibits a fi- 
nite superfiuid fraction due to quantum tunneling. How- 
ever, this superfluid fraction can be so small that the 
system can be considered virtually insulating. Although 
phase diagrams make rigorously sense in the thermody- 
namic limit, it should be taken into account that linear 
dimension of current experimental realizations of the sys- 
tem under concern is of the order of a few hundred sites. 
It is hence important to consider finite-size effects, which 
we demonstrate to be quite relevant especially in the par- 
tially compressible phases, and to depend significantly on 
the impurity density. 

The plan of the paper is as follows. In Sec. [Til we 
describe the model under investigation and introduce 
the superfluid fraction as an important parameter in 



the characterization of the phase diagram thereby. In 
Sec mil we recall the site-dependent Gutzwiller approach 
and provide analytic expression for the superfluid frac- 
tion and the flux induced by an infinitesimal velocity field 
in this framework. Section IIVI is devoted to the phase 
diagram of the system. First of all we discuss how the 
compressibility scenario argued in Ref. 22!| is captured by 
the site-dependent mean- field approach. In Sec. lIV Al we 
provide an analytical form for the boundaries of the fully 
incompressible insulating lobes. Moreover we discuss how 
the excitation spectrum of the system |5| is modified by 
the presence of the noninteger insulating phases charac- 
terizing of the Bose-Hubbard model with binary disorer. 
Section IIVBI discusses the superfluidity of the partially 
compressible phase, in relation to the quantum perco- 
lation phenomenon. Finally finite-size effects are inves- 
tigated in Sec. IIVCI for ID and 2D systems. This pa- 
per also contains a rich Appendix section where several 
interesting analytical results are provided for the site- 
dependent Gutzwiller approach. In particular in Sec. El 
we clarify the connection between the mean-field Hamil- 
tonian and the dynamical Gutzwiller equations. In Sec.lBl 
the analytic formula for the boundaries of the incom- 
pressible lobes is derived explicitly. Also we provide an 
useful inequality and clarify the connection between such 
a formula and similar results derived in single site mean- 
field approaches [H, [IH, [H [13, IMi . Finally, the superfluid 
fraction and the flux across neighbouring sites is derived 
analytically as a function of the mean-field order param- 
eters alone in Sec. [Cj A particularly simple expression 
applying in the ID case is also provided. 



II. THE MODEL 

The system under investigation is described by the 
Bose-Hubbard Hamiltonian 
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The on-site bosonic operators a^, a, and rij — a^aj de- 
stroy, create and count particles at lattice site j, respec- 
tively. The geometry of the A/-site lattice is described 
by the adjacency matrix A, whose generic element Aij 
equals 1 if sites i and j are nearest neighbors, and oth- 
erwise. The parameters U and T are the on-site repulsive 
strength and the hopping amplitude across neighboring 
sites, and, from the experimental point of view, they are 
related to the scattering length of the alkali atoms form- 
ing the bosonic gas and the strength of the optical lattice. 
We will be considering a binary random distribution for 
the local potential Vj, namely 



pi'"]) = PoS{vj - A) + (1 - Pq)S{v.j) 



(2) 



This choice is meant to account for the presence of 
A^imp = MpQ atoms of a second species trapped at ran- 
domly determined sites by a strong quench in the relevant 



hopping amplitude [a, [l^. The parameter A measures 
the strength of the interaction between these frozen im- 
purities and the bosons described by Hamiltonian ^ . In 
most of the following discussion we assume A < U. The 
general case can be worked out straightforwardly, and it 
is briefly discussed in Sec. IIVAI 

As it is well known pl, on a homogeneous lattice, 
Vj = 0, the zero-temperature phase diagram of the BH 
model described by Eq. ^ comprises an extended su- 
perfluid (SF) region and a series of Mott-insulator (MI) 
lobes. The SF phase is gapless, compressible and charac- 
terized by nonvanishing superfluid fraction. Conversely, 
the MI phase is gapped, incompressible, and character- 
ized by vanishing superfluid and condensate fractions. 
The presence of random potentials is expected to induce 
a further Bose-glass (BG) phase which, similar to MI is 
not superfluid, but, similar to SF, is gapless and com- 
pressible [1[. Recently it has been shown that in the 
case of uniformly box-distributed disorder such a phase 
can be captured by a multiple-site mean-field approach 
|13 . lis . Il4l , both on one- [IJ, [lj| and two-dimensional 
lattices [1J|. In the latter case the presence of the har- 
monic trapping potential typical of experimental systems 
was also taken into account. 

The superfluid fraction is estimated as the response of 
the system to an the infinitesimal velocity field imposed 
on the lattice. In the general case such a field is de- 



scribed by the antisymmetric matrix Bi 



-Bj i having 



nonzero elements only across neighboring sites. In the 
reference frame of the moving lattice the Hamiltonian of 
the system has the same form as Eq. ([T]) except that Aij 
is substituted by Aij exp{i9Bij), where is a scalar re- 
lated to the modulus of the velocity field [see e.g.[25l.[2a|. 
The superfluid fraction is often defined as the stiffness 
of the system under the pha se variation imposed by the 
velocity field [see e.g. [23, U^ 
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where Eg and Eq are the ground-state energies of the sys- 
tem when the lattice is moving and at rest, respectively, 
while N is the total number of bosons in the system. It 
should be noticed that Eq. ([3]) is properly a fraction^ i.e. 
a quantity with values in the interval [0, 1], only in simple 
situations, such as a homogeneous velocity field. More in 
general, the condition max(|i?,ij|) < 1 ensures that the 
superfiuid fraction does not exceed 1. The imposition of 
a velocity field induces a flux across neighboring sites i 
and j of the form 



J,, = ztAy(vI/|e-'^^-a]a, - e'«^-a|a,|vl/), (4) 



which clearly vanishes for 6 = 0. In the following we 
will show that /s = only if J^ij = for any pair of 
neighbouring sites. 



III. MEAN-FIELD APPROXIMATION 

The results we are going to illustrate are obtained in 
the widely used site-decoupling mean-field approximation 
[29l . [so] ■ Despite this approach cannot capture the cor- 
rect behavior of the spatial quantum correlations, it pro- 
vides a qualitatively satisfactory picture of the phases of 
strongly correlated systems, even in the presence of spa- 
tial inhomogeneities arising from the harmonic confine- 
ment typical of experiments |3ll| or from superimposed 
disordered potentials [£, jj, [lj| . 

In the strongly correlated regime the state of the sys- 
tem is expected to be well approximated by a Gutzwiller 
product state 
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where |51) is the vacuum state, aj\Q) = 0. A time- 
dependent variational principle similar to that illustrated 
in Ref. [321] results in a set of nonlinear dynamical equa- 
tions for the expansion coefficients Cj y [33| . It can be 
shown that finding the minimum-energy stationary state 
(fixed-point) of such equations is equivalent to finding 
the ground state of the mean-field Hamiltonian 



- ^(7^al+7^*a^). 
subject to the self-consistent condition 
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Such a mean-field Hamiltonian is usually derived by in- 
troducing the decoupling assumption 



aja. 



uj uj — a*aj in Hamiltonian ([T]) [30]. These issues will be 
briefly discussed in Appendix [X] The parameter /i ap- 
pearing in ([7]) is the so-called chemical potential, which 
comes about due to the fact that the mean-field Hamil- 
tonian ([6]) does not preserve the total number of bosons, 
unlike Eq. ^. 

The mean-field order parameters ai (with 6 = 0) allow 
the characterization of the quantum phases of the system, 
as far as compressibility is concerned. In particular, on 
a homogenous system the (site-independent) a^ is zero 
in the incompressible insulator and finite in compressible 
supcrfluid phase [30| . On inhomogcnous lattices a further 
situation can in principle occur, where ai ^ only on a 
fraction of the lattice sites. Following Ref. [24] we will 
deflne such situation as partially compressible, as opposed 
to the fully compressible and fully incompressible phases 
corresponding to the MI and SF. 



We observe that in this framework the one-body den- 
sity matrix has a rather simple expression, 



Pi 



(*]a|aj]*) _ 5,,{^j\n,\^Pj) + il- S^j)a*aj 



N 



N 



so that the condensate fraction, i.e. the largest eigenvalue 
of Pi j [3J] , can be estimated as 



f^'^N^ 



(10) 



Such a form shows that, at the mean field level, the con- 
densate fraction vanishes only for fully uncompressible 
MI phases. Actually fr has been used a convenient order 
parameter in Refs [l^, [l3, J_4, _35] . 

As one expects, the presence of the velocity field in- 
volved in the evaluation of the superfluid fraction mod- 
ifies the self-consistcntly determined mean-field param- 
eters defined in Eq. ([5|). It is easy to show that in 
the homogeneous case the superfluid fraction defined in 
Eqs. ^ equals the right hand side of Eq. ^U^. That is, 
the superfluid and condensate fraction coincide and can 
be equivalently employed for characterizing the phase di- 
agram of the system. In the general case a perturbative 
approach, carried out in Appendix [Cl shows that the su- 
perfluid fraction is 



1 



fs = ^l^Av'^>jiB^3-'i^^+^jy 



(11) 



where the real and positive a" are the mean-field order 
(7) parameters for 9 = 0, and the real phase factors (pj de- 



pend on the a*^ according to Eq. (|C7p . At the first per- 
turbative order in 9, the mean-field order parameter are 



aj = a'j eyip{i9(j)j), 
while the flux defined in Eq. fl]) becomes 
Jij- = 29ta°a'^jA,j[B,j - 4>, + > 



(12) 



(13) 



Clearly the superfiuid fraction in Eq. (|11|) vanishes only 
if each of the terms in the sum vanishes, i.e. if all of the 
fluxes in Eq. (|13p are zero. Expectedly, this happens in 
the MI phase, where a^ = at every site. The same 
can happen under more general conditions. Indeed, it is 
sufficient that (pi — 4>j = Bij whenever Q;°a° 7^ 0. This 
is precisely what happens in a BG phase, where fs — 
despite the system is compressible. 

It is interesting to note that on ID systems, owing to 
the conservation of fiux, the evaluations of Eqs. (fT5|) and 
(jlip does not require the determination of the phases (pj 
via Eq. (|C7[) . Indeed, as illustrated in Appendix [Ul one 
gets 



Jjj+l - -Jj+lj - J -2t9iY^ n n 



(14) 




FIG. 1: Expected phase diagram for a d dimensional lattice. 
The solid and dashed black curves are the (mean-field) bound- 
aries involved in Eq. (|16|l . These delimit three phases, as far as 
compressibility is concerned. The uncompressible lobes (dark 
gray), the partially compressible regions (light gray) and the 
fully compressible region (white). The uncompressible regions 
are labeled by the relevant filling. The partially compressible 
regions are labeled by the local potential of the favourable 
sublattice (see text for more details). The data points have 
been obtained as described in Sec. IIV Al for a very large ID 
lattice (M = 10^) and two impurity densities, po = 0.2 (white 
circles) and po — 0.6 (black circles). Note that both data sets 
agree very well with the expected density-independent ana- 
lytic resuh, Eq. (ITS)) . 
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(15) 
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It is clear from Eq. P^ that J^ = and /s = as soon 
as one of the local mean-field parameters vanishes. This 
explicitly shows that the partially compressible phase is 
insulating, and hence Bose-glass, in ID system, as al- 
ready mentioned [23, 123] • 



IV. PHASE DIAGRAM 

As discussed in Ref. [22|, the zero-temperature phase 
diagram of the system with binary disorder can be easily 
inferred from that of a homogeneous system, at least as 
far as compressibility is concerned. Indeed, independent 
of the impurity density Nimp/M, in the thermodynamic 
limit M — > oo a finite fraction of the disordered system 
consists of arbitrarily large regions of uniform local po- 
tential (Lifschitz rare regions, see [22, [23|). The bulk of 
these regions will behave as a homogeneous lattice, un- 
dergoing a transition at the analytically known critical 
value 



U 2d 
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g(x)=^"'LxJ)(N-x) ^^g^ 



where [a;J denotes the largest integer smaller than x, 
\x~\ — lx\ + 1, d is the lattice dimension and v is the 
local potential within the homogeneous regions, which 
can attain the values and A. It is clear that the sys- 
tem is fully compressible or fully incompressible only if all 
of the above region exhibit the relevant property. This 
means that the region of the phase diagram 
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is fully compressible, while the incompressible lobes cor- 
respond to the region 



^<B. 
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Fig. [T] shows the phase diagram of the system for A = 
0.5 U. The fully compressible and fully incompressible 
regions correspond to white and dark gray shading. As 
in the homogeneous case, the fully incompressible lobes 
correspond to plateaus of the system filling. Interestingly, 
the binary disorder causes the appearance of non-integer 
critical fillings of the form [/x/f/] —po, where po is the im- 
purity density [2l|, [22, [23| . The incompressible MI phases 
will be discussed in detail in Section llV Al The light gray 
region of Fig.[Tl enclosed between the boundaries Bi and 
B2, is of course compressible, yet it contains arbitrarily 
large incompressible regions. Therefore it is referred to 
as partially compressible [22] . 

This apparently simple scenario requires some clarifi- 
cations when the possible superfluidity of the system is 
concerned. As it can be understood from Eq. P^ . the 
system can sustain a superfluid flow only along a path 
where the mean field parameters aj are not vanishing. 
Hence, as expected, the system is not superfluid in the 
incompressible regions of the phase diagram. Conversely, 
in the fully compressible region the mean-field parameters 
are nonzero almost everywhere, and a finite superfluid 
fraction is expected. The most interesting regions are 
the partially compressible ones, which can behave as BG 
phases, as it will be discussed in Sections IIV BI and IIV CI 



A. Incompressible phases 

The boundaries of the incompressible Mott lobes have 
been derived in Ref. [2l[ based on a single-site mean-fleld 
effective theory. Reference [23] reports more quantita- 
tive results ensuing from strong-coupling perturbative ex- 
pansions, which are further supported by densit y rn atrix 
renormalization group simulations. Reference [23| also 
provides results based on strong-coupling perturbative 
expansions, as well as on exact diagonalization of small 
ID systems. Here we discuss site-independent mean-field 
results and show that, unlike effective single-site mean- 
field theories, they do not depend on the impurity density 
Pq, as it is expected. 



We first of all observe that in the so-called atomic limit 
t ^ the ground-state of Hamiltonian ([1]) is a product 
of on-site Fock states. That is, Eq. ^ applies exactly 



spond to different impurity densities, TVii 



0.2M and 



with 
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where Vj = maxjO, [/i — Vj~\} and the 



Mmp = 0.6M. Both data sets show a very good agree- 
ment with the analytic result in Eq. psp and, as ex- 
pected [22|, |23|, exhibit no dependence on the impurity 



chemical potential /i is determined by the constraint on density po = Ni^p/M. As briefly recalled in Appendix iBl 
the total number, N — J2i '^j- Recalling that the local 



potential Vj is A at iVimp = Po M randomly placed lattice 
sites and at the remaining AI — A^imp sites, it is easy 
to conclude that the total number of bosons is zero for 
— oo < /i < 0, and subsequently grows stepwise with 
increasing chemical potential. The staircase function is 
easily determined if < A < [/. In this case the rises of 
the steps occur at /ifc(x) — k + x, with k = 0, 1, . . . , cxd 
and X ~ 0, A. The height of the steps, i.e. the total 



population, is A^ = M{k + 1) 



Nu 



for fj.k{0) < /z < 



/ifc(A) A^ = M{k + 1) for ^fe(A) < fj. < Mfc+i(0). In the 
first case the wavefunction ([5]) is such that \4>j) = |A:) at 

= |fc -I- 1) at the 
'kl is the local k- 
- |A: + 1) at every 



the A^inip sites with Vj — A, and |i/'j) 
remaining sites, where \k) = {a-)'^\^)/ 
th Fock state. In the second case \^pj) 
lattice site. 

It is straightforward to check that these states diago- 
nalize the mean-field Hamiltonian ([6]) subject to the self- 
consistency constraint ^ also for any i > 0, although 
they do not always represent the mean-field ground-state 
of the system. As it is illustrated in Sec. [BJ this is true 
only within the regions of the n/U-t/U phase plane de- 
scribed by < t/U < |Amax|~^, where Amax is the maxi- 
mal eigenvalue of the matrix 



A = DA. 
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A is the adjacency matrix of the lattice and the func- 
tion B appearing in the diagonal matrix D is defined in 
Eq. P^ . The numerical diagonalization of A at different 
values of /i shows that the above discussed plateaus ex- 
tend over lobe-like regions with alternatively noninteger 
and integer fillings. 

Since in both cases aj = we classify these phases 
as incompressible Mott insulators. As we have discussed 
above, the superfluid fraction cxpcctedly vanishes. How- 
ever, integer- and fractional-filling insulating phases are 
distinguished by the correlation with the underlying lo- 
cal potential. The former are homogeneous despite the 
presence of such potential. The latter are clearly charac- 
terized by a disorder directly related to that in the loca- 
tion of the impurities described by the local potential Vj . 
We observe that these two uncompressible phase are ex- 
pected to exhibit different excitation spectra, which are 
relevant experimental quantities [5]. More to the point, 
the excitation spectrum of the integer-filling Mott phases 
will be characterized by three peaks at U — A, U and 
U + A. As to the noninteger-fiUing lobes, the peaks are 
expected at A, U and, for fillings larger than 1, at 2U—A. 

The data points in Fig. [1] have been obtained by eval- 
uating the maximal eigenvalue of the matrix A as a 
function of the chemical potential for a ID lattice com- 
prising M = 10^ sites. Black and white circles corre- 



imp/ 

the same result can be equivalently obtained by studying 
the matrix F = yDAyD which, unlike A, is symmetric, 
F = F*. 

We emphasize that mean-field results based on effec- 
tive single-site results can be recovered by averaging the 
above matrices over disorder. We first of all observe that 
the disorder-averaged version of A results in the critical 
boundaries 
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i.e. precisely the same result obtained by the so-called 
sim,ple man's mean-field theory of Ref. [2l|. Note that 
Eq. (pO)) is indeed a simplified result, in that it depends on 
the impurity density poi contrary to expectations [22ll23|. 
We also observe that equivalent simplified approaches 
have been adopted in earlier papers [ll|, |l3| , albeit with a 
different disorder distribution, and date back to the sem- 
inal work by Fisher et al. [ll] ■ A perhaps more structured 
effective single-site theory is obtained from the disorder 
averaged version of the matrix F matrix defined above. 
Indeed the boundary of the aj = phase ensuing from 
such matrix in the case of uniformly distributed disorder, 
p{vj) = Q{vj + A/2)0(A/2 — Vj), is very similar to that 
provided by the stochastic m,ean-field theory described in 
Ref. [lal- However, it is quite clear that in the case of 
the binary distribution in Eq. ^ such boundary again 
depends on the impurity density po, just like in Eq. (^0)1 . 

In the above general discussion we assume < A < U. 
For A > U the arrangement of the first few lobes changes 
straightforwardly. For instance, for U < A < 2U the 
unitary-filling lobe disappears, the basis of the lobe at 
filling 1 — Po extends in the whole interval < /i < C/ 
and the interval U < /.i < A provides the basis for a 
lobe with filling 2 — 2po. Note that if po = 1/2 this 
last lobe has unitary filling, although this results from 
averaging sites at filling and 2. Hence the unitary- 
filling incompressible phase changes from homogeneous 
to disordered as A becomes larger than U. This was 
observed in an early work, where however the disordered, 
incompressible insulating phase is identified as a Bose- 
glass phase 20] . 

Clearly, the phase diagram of the homogeneous lattice 
is recovered for A = 0, while the case A < can be 
mapped on the repulsive case with the suitable number 
of impurities, A -^ |A|, A^imp ^ M — A^imp- 



B. Partially compressible phase, Bose-glass 

In order to discuss the situation in the partially com- 
pressible regions it proves convenient to introduce the 
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FIG. 2: Power law decay of the superfluid fraction as a func- 
tion of t/U . The left panel refers to a ID lattice comprising 
M = 400 sites. The different data sets have been obtained by 
varying the impurity density between 0.4 and 0.7. The right 
panel corresponds to a 2D system, also comprising M = 400 
sites, where the impurity density varies between 0.005 and 0.1. 
In every case the data sets are well described by a straight 
line having an integer slope which is found to coincide with 
the percolation length L (see Fig. [31). In the ID case (left) we 
observe lengths from 3 to 14. In the 2D case L = 3, 4, 5, 6. 







FIG. 3: A sketch illustrating the notion of percolation length. 
A 30 by 30 square lattice contains 14 favourable sites, sig- 
nalled by the black squares. Due to their very low density 
they do not form a cluster spanning the lattice. The dashed 
line signals the shortest path allowing to "wade" through the 
lattice by stepping on favourable sites. The longest jump has 
to be taken between the first and the second favourable site 
from the left. The gray shading demonstrates the effective 
increase in the density of favourable sites caused by the long- 
range connectivity introduced by quantum tunnelling. The 
range of such connectivity is 4 (left panel) or 5 (right panel) 
lattice constants. The shading in the rightmost panel shows 
that the longest jump in the dashed path measures 10 lattice 
constants. Hence L = 10 for this favourable sublattice. 



notions of favourable and unfavourable sublattices for 
bosons added to the system, as determined by the com- 
petition by the local boson-boson interaction strength U 
and the local impurity potential A < U. Clearly, in the 
absence of bosons the impurity-free sublattice, Vj — 0, is 
energetically favourable. If some bosons are introduced 
in the system, they will prefer the impurity- free sites, as 
far as the local energy is concerned. However, at some 
fillings, the presence of bosons in the impurity- free sub- 
lattice could make the impurity sites more convenient 
energetically. Thus the favourable sublattice coincides 
with the impurity-free sublattice or with the impurity 
sublattice depending on the filling or, equivalently, on 
the chemical potential. Typically, the impurity sites are 
preferred in the partially uncompressible regions to the 
left of the integer filling uncompressible lobes (marked by 
A in Fig. [1]) , whereas the impurity- free sites are preferred 
in the remaining regions (marked by in Fig.[T]). We call 
unfavourable the sites of the lattice not belonging to the 
favourable sublattice. 

The partial compressibility of the regions we are con- 
sidering arises from the presence of arbitrarily large un- 
favourable regions behaving like homogeneous systems. 
Hence there are arbitrarily large sublattices hosting an 
uncompressible phase characterized by aj =0. However, 
owing to the hopping term, not all of the unfavourable 
sites are characterized by vanishing order parameter. 
This allows for the formation of a cluster of aj > 
sites spanning the entire lattice, and therefore capable 
of supporting a superfluid flow, even if the favourable 
sites do not percolate throughout the lattice. That is to 
say, the superfluidity is not related to the percolation of 
the favourable sublattice, but rather to the percolation 
of the sites where a^ > 0. The latter is made possible by 
the quantum tunneling effect, which allows for the bridg- 



ing of possibly disconnected clusters of favourable sites. 
We have studied this phenomenon at fi — A — 0.5, i.e. 
where the partially uncompressible regions of the phase 
diagram extend down to vanishing hopping amplitude. 
As it is shown in Fig. [51 both in ID and 2D we observe 
a behaviour of the form /g oc (t/t/)^, where L is the (in- 
teger) percolation length, i.e. the length of the longest 
"bridge" among those necessary to turn a disjoint im- 
purity distribution into the shortest cluster spanning the 
whole lattice (see Fig. [3]) . Note that the same po can give 
a different L due to finite-size fluctuations. This is espe- 
cially clear on ID system, where one expects L = oo in 
the thermodynamic limit, independent of po- However, 
once an L is determined by the actual finite-size realiza- 
tion of the disordered system, it dictates the behaviour 
of /s as discussed above and demonstrated in Fig. [3 

Now on one-dimensional lattices the percolation length 
corresponds to the size of the largest homogeneous clus- 
ter in the unfavourable lattice, which becomes arbitrar- 
ily large in the thermodynamic limit, independent of the 
impurity density. Hence in ID the partially compressible 
phase is expected to be insulatin g fo r any finite impurity 
density as observed in Refs. |2a.l23j. 

Conversely, on higher dimensional lattices the super- 
fluid fraction is finite at any impurity density, although 
it can become extremely small. This is due to the fact 
that, unlike the one-dimensional case, L never diverges in 
the thermodynamic limit. This can be understood by ob- 
serving that quantum tunneling introduces a long-range 
connectivity between the possibly disjoint clusters form- 
ing the favourable sublattice, which effectively increases 
the density of favourable sites. It is then clear that a suf- 
ficiently large hopping amplitude can bring the effective 
density of favourable sites above the (finite) percolation 
threshold, so that an effective spanning cluster is formed 



[301 . This concept is illustrated by the light gray regions 
in Fig. El 

The above discussion is valid in the thermodynamic 
limit. As we illustrate in the following, strong deviations 
from the expected behavior can be observed on finite-size 
lattices, even for fairly large sizes. In particular, on ID 
systems a finite superfluid fraction can be observed in 
the partially compressible regions of the phase diagram. 
Conversely, on higher dimensional system, the expectedly 
finite superfluid fraction in these region may become so 
small that the phase can be considered insulating for any 
practical purpose. This agrees with the intuitive notion 
of a glass as an extremely viscous fluid. 



C. Finite-size effects 

Of course a phase diagram is rigorously defined only in 
the thermodynamic limit, M — > oo. However, the exper- 
imental realizations of the Bose-Hubbard model, based 
on ultracold atoms trapped in optical lattices, are far 
from such a limit, especially in the presence of disorder. 
Indeed, the occurrence of arbitrarily large regions with 
uniform local potential becomes extremely improbable, if 
not impossible. The higher the lattice dimension d, the 
more serious this problem. Moreover, as demonstrated 
by Fig. [21 the superfluid fraction can be so small that 
the system can be considered virtually non superfluid. 

In order to demonstrate the relevance of these effects 
we have carried out numerical simulations for ID and 2D 
lattices comprising A/ — 961 sites. In both cases we have 
adopted periodic boundary conditions, and a constant 
velocity field parallel to a coordinate direction. The re- 
sulting phase diagrams, where we took into account that 
vanishingly small superfluid fractions can be considered 
zero, are shown in Figs. [H and [5l The first thing to be no- 
ticed is that finite-size effects do not dramatically affect 
the boundaries of the fully incompressible Mott lobes, 
especially in one dimension. A slight dependence on the 
impurity density po is observed, at variance with the ther- 
modynamic limit result. In general, as it is explained in 
Appendix [Bl the finite size lobes "enclose" those in the 
thermodynamic limit. Significant finite size effects are 
instead evident in the partially compressible phase. As 
we mention above, in the thermodynamic limit one ex- 
pects this phase to be insulating in ID, and superfluid for 
d> 1. Actually, as it is clear from Figs. [H and [5l these are 
the predominant characters of the partially compressible 
phases also on flnite-size lattices. However, on ID lattices 
signiflcant portions of the partially compressible phase 
are superfluid. Conversely, on 2D lattices small partially 
compressible regions surrounding the Mott lobes exhibit 
exponentially small superfluid fractions, so that they can 
be considered virtually insulating. While the expected 
behaviour in the thermodynamic limit is expected to be 
independent of the impurity density po, this parameter 
strongly affects the finite size deviations. In particular, 
for small impurity densities the partially compressible re- 




FIG. 4: Phase diagram as determined from a numerical sim- 
ulation on a ID lattice comprising M — 961 sites and con- 
taining Nin-ip = 200 (left) and Mmp ~ 600 (right) randomly 
located impurities. In both cases A = 0.5 U, as in Fig.[T] The 
density plot represents the superfluid fraction as specified by 
the colorbar. The dark grey areas enclose the region where /s 
is smaller than 1% of its largest value in the entire examined 
area. The light grey areas are the fully incompressible Mott 
lobes determined as described in Sec. IIV Al The solid and 
dashed black curves are the boundaries involved in Eq. 1161 
They are the same as in Fig. [T] 
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FIG. 5: Phase diagram as determined from a numerical sim- 
ulation on a 2D lattice comprising M = 31x31=961 sites 
and containing Nimp = 200 (left) and Nimp = 600 (right) 
randomly located impurities. In both cases A = 0.5 U, as in 
Figs. [Hand H 



gions marked A in Fig. [I] tend to be more insulating than 
thosed marked 0. The converse occurs at large impurity 
densities. 

We conclude by emphasizing that the favourable sub- 
lattice does not percolate at neither of the chosen impu- 
rity densities. Hence, the fact that /s > in the partially 
compressible phases is to be attributed to quantum tun- 
neling effects. 



V. SUMMARY 

In this paper we address the phase diagram of the 
Bose-Hubbard model describing ultracold bosonic atoms 
loaded in an optical lattice containing static random lo- 
cal impurities. These are fermionic atoms whose hopping 
amplitude has been quenched to extremely low values. 
We employ a site dependent Gutzwiller scheme, analyz- 
ing both ID and 2D lattices. On the one hand we show 
that this approach confirms that the phase diagram of 



the system does not depend on the density of impurities 
and can be easily derived from the phase boundary of 
the homogeneous case, at least with respect to the com- 
pressibility of the system. We show that the boundaries 
of the insulating fully incompressible region of the phase 
diagram are strictly related to the spectral radius of two 
(block) tridiagonal matrices. We discuss the expected 
modifications in the structure of the experimental excita- 
tion spectrum of the system occurring due to the presence 
of the noninteger-fiUing incompressible lobes appearing 
in the phase diagram in the presence of binary disorder. 
Also we provide exact formulas for the superfluid fraction 
and the fluxes induced by an infinitesimal velocity field, 
showing that these quantities ultimately depend on the 
mean-field order parameters alone. These formulas allow 
us to investigate the superfluid nature of the partially 
compressible regions appearing in the phase diagram of 
the system owing to the binary disorder. We discuss 
the experimentally relevant finite size effects, showing 
that they strongly depend on the impurity density and 
mainly affect the boundaries of the partially compress- 
ible regions. We show that on one-dimensional system 
quantum percolation causes the appearance of superfiuid 
domains within these regions. On the other hand, on 
higher dimensional lattice the in-principle superfluid par- 
tially compressible regions contain domains that can be 
considered virtually insulating, due to the extreme small- 
ness of the superfiuid fraction. 

Last but not least, in the appendices we provide the 
explicit derivation of the analytic results employed in the 
paper, whose validity is not limited to the case of binary 
distributed disorder. In particular, we discuss the rela- 
tion the site-dependent Gutzwiller approach and effective 
single site mean-field theories. 
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APPENDIX A: NORMAL MODES OF THE 
GUTZWILLER DYNAMICS 



As we mention in Sec. IIIIl a variational principle anal- 
ogous to that described in Ref. [30] results in the set of 
dynamical equations for the expansion coefficients of the 



Gutzwiller state ([5]) [^ 

U 



ic 



JV - -^V{y- \)Cjy+V^VCjy 



- t {j*\Jv \- Icj y+i 4- ^j\fvcj j,_i) (Al) 

OO 

7j = ^Aj/jCk/,., g/i =^ y^ \/t^ + 1 c* ^Cj ^+1 (A2) 

h 1^=0 

It is straightforward to check that the norm of each on- 
site Gutzwiller factor, and the (average) total number of 
bosons in the system 



(V',i^.> = Ei 



^ = EE 



(A3) 



are conserved by the dynamics. 

The fixed-point condition Cj ^ = Q for Eq. (|A1[) results 
in the set of equations 



= 



-(■J + i^^i^ - 1) + ("j - /^) 



L-J 1/ 



- t {%j\/v+ Icj y+i + -/jv^Cj ^„i) (A4) 

where /i and {ej} are M-l-1 Lagrange multipliers ensuring 
that the total number of bosons and the norms of the 
Gutzwiller factors |'0j) equal the desired values. 

But Eq. (|A4[) is nothing but the eigenvalue equation 
for the on-site mean-field Hamiltonian ([7]) projected onto 
the generic Fock state at site j, Iv) = {aA'^\Q)/ViA, 



Mn, 



iHj 







(A5) 



where we recall that ah\^) — for all /I's. Note in- 
deed that Eq. (|A2|1 . which must be considered part of 
Eq. lAll is exactly equivalent to the self-consistency con- 
straint specified by Eq. ([8]). By comparing Eqs. (|A4p 
and (|Aip we see that the solutions of the mean-field equa- 
tions are normal modes of the Gutzwiller dynamics such 
that Cjy{t) = e-'-^^^^+f'^hji^iO). The fixed-point nature 
of these solutions becomes clear when one considers the 
relevant expectation values on Hcrmitian operators, cor- 
responding to observable quantities. It is easy to verify 
that number-conserving Hcrmitian operators, such as for 
instance aj^ah + aj^ak, produce time-independent expec- 
tation values, whereas non number-conserving Hcrmitian 
operators, such as aj^ -\- ah, produce expectation values 
oscillating around 0. Note that the latter should be iden- 
tically zero, since the original Hamiltonian ([T]) — unlike 
its mean-field counterpart — commutes with the total 
number of bosons. This result is recovered after time- 
averaging the expectation values. 

Note finally that E = (*|i?|*> = (*|W + M^|*), 
and that Ti can be obtained from H — fj,N by assum- 

a*ak for j ^ fc, where 



a]afc 



(V'ikilV'i 



ttj-afc 



ing that a}:ah 

a, = (*|aj|*) 

As it is shown in Refs. [3j, |37| the TDVP approach 
based on Glauber's and SU(M) coherent states instead 



of those in ([5]) gives the discretized Gross-Pitaevskii equa- 
tions for the Bose-Hubbard model H]). Interestingly, the 
Gutzwiller mean-field states in Eq. in ([5]) reduce to 
Glauber's coherent states for C/ ^ 0, as it is clear from 
the form of the mean- field Hamiltonian ([6]) [32] . 



APPENDIX B: MOTT PHASE BOUNDARY 



Restoring the site label, recalling the definition of 7„ 
Eq. (HI), I/, Eq. ^^ and B, Eq. ^ one gets 



t K?— 1 I f^ ^"^ 



u \ u 



U 



771 m' \*-^rn' 



/ ^ -^'m m' \P"n 



(B7) 



In this appendix we show that the critical boundary of 
the (mean-field) Mott phase, [ctj — for every j) is the 
inverse of the maximal eigenvalue of the matrix A defined 
in Eq. ^. 

Since aj = everywhere, we also have 7^ = at every 
site, according to Eq. ([5]). Hence the mean- field Hamilto- 
nian (l6|) is the sum of the on-site Hamiltonians in Eq. ^ 
and the the ground state ([5]) is bound to be a product of 
local Fock states 



(at) 



t^" 



■|0), 



The relevant on-site energy is 

1) 



^ ( 



' u ' 



li)y 



(Bl) 



(B2) 



Hence (^|aj |^) = at every site, and the self-consistency 
constraint ([5]) is satisfied. Note that Eq. ([5]) defines a 
map, since any set of (possibly nonzero) aj determine 
a set of local ground states |-0j) via Hamiltonian ([6]), 
which in turn determine a new set of Uj . This is by def- 
inition a fixed point of the map when it coincides with 
the original set, which is exactly what happens for the 
configuration under examination, aj — Q, for any choice 
of the Hamiltonian parameters. However, the stability 
of such "trivial fixed point" does depend on the Hamil- 
tonian parameters. Specifically, the fixed point is stable 
only if the maximal eigenvalue of the linarized version 
of the map is smaller than 1. In order to linearize the 
map we assume that \aj\ <C I, which yields \^j\ <^ 1 and 
treat the (mean-field) kinetic term in Hamiltonian ^ as 
perturbative. Dropping for a while the site label we get, 
up to the first perturbative order. 



IV^) - I^W) + I^A^i)) 



|^W)^-i 



E 



{v'Y^a^ + 7*^1^) I 



(B3) 



(B4) 



\^) = \v)--^ \v + \) ^^^1^.-1) (B5) 

where \v)t v and e^ are defined in Eqs. (|BI|) and (|B2p . 
Hence 



'\^\aW) = -^7 



i^ + l 



-7 






(B6) 



where A is the same as in Eq. ([T9|) . Recalling the criteria 
for the stability of linear maps, the fixed point (am) = 
(equivalent to 7™ = 0) is stable whenever 



t 
— < 



1 



(B8) 



where A^ax is the eigenvalue of A with the largest mag- 
nitude. 

We note that despite A* 7^ A, the spectrum of this ma- 
trix is real. This can be explained by observing that the 
eigenvalue problem Ax = Ax is equivalent to Fy — Ay, 
where yj — D~^Xj and F — ^/DA^/D = F*. Recalling 
that the maximal eigenvalue of a matrix coincides with 
its 2-norm one can derive a lower bound for the critical 
hopping to interaction ratio. Indeed 



I 



I 



> 



I 



2d min B 



l^l|2P||2 

u 



(B9) 



Note that in the case of binary distributed disorder 
Eq. (jB9[) coincides with Eq. p^ . Hence, the finite-size 
lobes always enclose their thermodynamic counterparts. 
It is interesting to observe that the above approach 
naturally suggests two disorder-averaged effective theo- 
ries. Indeed, one could trade the matrix elements of A or 
F with their averages over disorder. This would give 



t 




1 





<- 





u 




2d 



dvp{v) B 



fJ- 



U 



for A and 



U 



(BIO) 



(BIl) 



for F. 

As mentioned in Ref. [I4I, Eq. (|B10|) gives the 
zero-temperature (analytical) phase diagram derived in 
Refs. [l|, [ill, [l7| in the case of v uniformly distributed in 
[—A, A]. Interestingly, Eq. (jBf ip gives a different result, 
very similar to that obtained by the stochastic mean-field 
theory recently described in Ref. [I5]. 

The integrations in Eqs. (|B10[) and (jBlip can be easily 
carried out analytically in the case of a binary distributed 
disorder, Eq. ([2]). In particular, it is quite straightforward 
to show that Eq. (JBIOP is equivalent to Eq. ([20|l . which 
does not exhibit the expected independence on the im- 
puriy density po = Ni^p/M in the thermodynamic limit, 



10 



as discussed in Sec. IIVI The same problem affects the 
boundary derived from Eq. (|B11[) . These results seem 
to suggest that — at least in the case of binary disorder 
— single-site mean-field theories are not able to capture 
the thermodynamic limit. 



APPENDIX C: SUPERFLUID FRACTION AND 
MEAN-FIELD APPROACH 

In this section we derive the analytic expression for 
the phases (pi introduced in Section Hill for the mean field 
study of the currents present in the system. Also, we 
obtain the analytic mean-field expression (|11[) for the su- 
pcrfluid fraction ([3]). 

We first of all observe that expanding the 6'-dependcnt 
terms in Eq. ([1]) one finds that even and odd contribu- 
tions are purely real and imaginary, respectively |28j . As 
a result, this is true also of the ground-state perturba- 
tive expansion, while only even terms contribute to the 
ground-state energy. 

In the Gutzwiller approximation, this implies that each 
factor of the perturbative expansion in Eq. ([S]) has the 
same alternating form: the even contributions are real 
and the odd contributions are imaginary. Hence, the 
same property holds also for the mean-field parameters 
am defined in Eq. ([8]). 

Therefore, at first order in 9 we can write 



aj = Oj eyip{i9(j)j) 



(CI) 



where a'j = (V'?|ail'0?) is the local mean-field parameter 
for 9 = 0, and (f>j S K. Plugging last result in Hamil- 
tonian ([7]) and keeping only first-order contributions one 
gets 

^j = -^i^j - 1)% + (vj - IJ-)nj - ta+Tj exp{idj) + c.c. 

(C2) 



where we denote 



^,= 



E.^, 



J«"i 



E^. 



a. 



(C3) 



Let us define hj = Oj exp(— z'iJj). The new creation and 
distractions operators aj and fij satisfy the same alge- 



bra of aj and aj and moreover luj 



Tij = ajUj. By intro- 
ducing these new operators the eigenvalue and the self- 
consistency equations become: 



EA^'j 



— {uj - l)nj + {vj - fi)nj 



<^] 



oE^. 



ji^j 



\^3 



exp(ii9j)('0j|dj|7/'j) = a'^ exp{i9(j)j 



(C4) 



(C5) 



The solutions to equations (|C4p and (|C5p can be directly 
obtained from the solutions of order in 9. In particular 
if iV'j) = En Cj,n(a^)"|0) is the on site wavefunction for 
6* = 0, to the first order in 9 the on site wavefunctions 



iv^°) = Ec^-^"(^.^)"io) = E^^-^"«""''^«)"io) (C6) 



with "di 



Moreover, as expected, there are no first- 



order contributions to the energy. Equation dj = 00j 
entails that 



Z^ 



Sj.'^Ajkal -Ajia° 



= Y.B,.al (C7) 



Equations (|C7|) provide an expression for the phases (pi, 
which can be plugged into equation (J13p . to obtain the 
currents to the first order in 9. 

Let us consider the the second order perturbation in 9. 
The contribution to the parameters aj can be written as 
aj = (a'y -|- 9'^S^j) e:x.p(i9(f>j). Within such approximation 
the on site Hamiltonian is 

Tij = l/2{nj — l)nj + {vj — iJ,)nj — ta'Jr', exp{i9(j)j) + c.c. 

(C8) 
with 

i \ i 



Equation (jC9p has been obtained by expanding in 9 and 
exploiting expression (jC7[) . Therefore we have to solve 
the self consistency equations: 



and 



where 



n''^+9'V^\ij,)=E,\tP,) (CIO) 



(^,|a,|V,) = (a° + 02^,) (Cll) 



-t{a+ + aj)J2Ajia° 



V, - -tih+ + h,) 



^-J2A,,a^{B,, + cP,-(b,)A (C12) 



Equations (JCIOI IC11|) is solved by a first order per- 
tubation theory in 9^. In particular denoting \tpj) = 
l^j) + ^^l^i) ^^^ Ej = E^ + 9^E^ (with |Vi°) given by 
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(|C6[) and Ef eigenvalue of the solution for 9 = 0), we 
have 



E^ = {i.^\V,\i.^), 



|^|) = -(i/5-i?^)-V,|^°) 



and 



(C13) 

(CM) 

(C15) 

(C16) 
withX,=t{i^]\{&,+&+){H°j-E^)-H~a, + ~a+)\i^']). The 
quantities ^j are obtained from Eq. (IC16P and Ej is 



inserting Eq. (jC14p into (JCISP we obtain 



^f 






(C17) 

We recall that the mean field Hamiltonian ([6]) is com- 
posed by the sum of the on site terms Hi and of a diag- 
onal term, which has been so far neglected, since it does 
not provide any contribution to the wave functions, and 
hence to the relevant expectation values. However, such 
term is not negligible in the evaluation of the system en- 
ergy. In particular, it provides a contribution of ^ E'l 
with Ef = a*t/2j2iAjiO:iexp{i9Bji) + c.c. For small 9 



we get Ef 



Ef 



9'^Ef with 



\ i i 

+ic,E^j»«? 

i 

from definition ([3| it is clear that 



(C18) 



(C19) 



It is interesting to observe that on ID systems there is 
no need to evaluate the phases 0j , the superfiuid fraction 
being determined by the a° alone. This can be proved 
by observing that on a ID lattice the current in Eq. (fT^)) 
is bound to have the same value across any couple of 
neighbouring sites j and j + 1. Recalling that Aij = 



3jj + l - 



'Si,j^i and Bi 



Oi,j + l ■ 



-Jj+ij ^J = -20ta5aO+i[0j 



5ij_i, one gets Jj.j+i = 
- 0J+1 — 1]. This yields 



2Mte 



and 



M ^ 



(C20) 



Pj - CPj+l 



-1 = - 



J 



29ta"^a"^^^ 



M ( "' I \ 



Plugging this result into Eq. (jC19|) gives 



fs = ^r 



M^ (^ 1 



J 



N l^^a°a«+i 



2t9N/M ■ 



(C22) 



We conclude by observing that the very same deriva- 
tion of Eq. (jC19[) can be carried out in the case of dis- 
crete Gross-Pitaevskii equations, which, as described in 
Sec. El ensue from assuming that the states in Eq. ^ 
are Glauber's or SU(Af) coherent states. The resulting 
equations have exactly the same form as Eqs. (|C7p and 
(|C19p . where the mean- field order parameter aj is sub- 
stituted by the corresponding coherent-state label Zj. 
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